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1 Introduction 

One hundred years ago, Albert Einstein completed his gen¬ 
eral theory of relativity, radically revolutionizing our under¬ 
standing of gravity. In this theory, gravity is no longer a 
force, but instead an effect of spacetime curvature. The mass 
and energy content of spacetime creates curvature that in 
turn dictates the behavior of objects in spacetime. Generally 
speaking, when a massive object is accelerating, it produces 
curvature perturbations that propagate at the speed of light. 
Such ripples in the fabric of spacetime are called gravitational 
waves (GWs). 

Forty years ago, Hulse and Taylor discovered the binary- 
pulsar system PSR B1913H-16 H]. Subsequent observations 
in the following years showed that its orbital period was grad¬ 
ually decreasing, at a rate that was entirely consistent with 
that predicted by general relativity as a result of gravitational 
radiation ijZjO. This provided the most convincing evidence 
of the existence of GWs and earned Hulse and Taylor the No¬ 
bel prize in Physics in 1993. 


*Con'esponding author (Xingjiang ZHU, email: xingjiang.zhu@uwa.edu.au; Linqing 
WEN, email: linqing.wen@uwa.edu.au) 


It was realized in late 1970s that precision timing obser¬ 
vations of pulsars can be used to detect very low frequency 
GWs (~ 10-®-10-^ Hz; Bllsl l. Based on a calculation by Es- 
tabrook & Wahlquist ||6l for the GW detection using Doppler 
spacecraft tracking, for which the principle is similar to pul¬ 
sar timing, Detweiler explicitly showed that “a gravitational 
wave incident upon either a pulsar or the Earth changes 
the measured frequency and appears then as a anomalous 
residual in the pulse arrival time" ’ l5l. In 1983, Hellings & 
Downs 13 used timing data from four pulsars to constrain 
the energy density of any stochastic background to be < 10“^ 
times the critical cosmological density at < 10“^ Hz (see, 
refs, mi for similar results). Almost simultaneously to these 
works, Backer et al. discovered the first millisecond pulsar 
PSR B1937-1-21 IfTOl . Because of its far better rotational sta¬ 
bility than any previously known pulsars and relatively nar¬ 
rower pulses, timing observations in the following years im¬ 
proved the limit on stochastic backgrounds very quickly and 
by orders of magnitude lfTTHT4ll . 

Measurements with a single pulsar can not make definite 
detections of GWs whose effects may be indistinguishable 
from other noise processes such as irregular spinning of the 
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star itself. By continued timing of an array of millisecond 
pulsars, i.e., by constructing a pulsar timing array (PTA), 
GWs can be searched for as correlated signals in the timing 
data. Hellings & Downs El first did such a correlation analy¬ 
sis to put limits on stochastic backgrounds. Romani M and 
Foster & Backer ED explored the greater scientific potential 
of PTA experiments: (1) searching for GWs; (2) providing a 
time standard for long time scales; and (3) detecting errors in 
the solar system ephemerides. 

There are three major PTAs currently operating: the Parkes 
PTA (PPt4]; ifTTlfTSl l was set up in 2004 using the Parkes ra¬ 
dio telescope in Australia; the European PTA (EPTaH IfTOl ) 
was initiated in 2004/2005 and makes use of telescopes in 
Erance, Germany, Italy, the Netherlands and the UK; and 
the North American Nanohertz Observatory for Gravitational 
Waves (NANOGra\ll JSOl) was formed in 2007 and car¬ 
ries out observations with the Arecibo and Green Bank tele¬ 
scopes. We summarize in Table [T]information about the tele¬ 
scopes used and the number of pulsars monitored by these 
PTAs. Recently the three PTA collaborations were combined 
to form the International PTA (IPTA0; HIBO l. Looking 
into the future, GW detection with pulsar timing observa¬ 
tions is one of major science goals for some powerful future 
radio telescopes, such as the planned Xinjiang Qitai 110 m 
radio telescope (QTT; which is a fully-steerable single-dish 
telescope EH) the Eive-hundred-meter Aperture Spher¬ 
ical Telescope (EAST; which is expected to come online in 
2016 |l25]|26l) in China, and the planned Square Kilometer 
Array (SKA) and its pathfinders (see ref. Il27l and references 
therein). Eigure [1] shows the distribution of IPTA pulsars on 
the sky along with those presently known millisecond pulsars 
that may be useful for PTA research with future telescopes. 


Table 1 Information about pulsar timing array projects. 



Telescope 

Diameter (m) 

Country 

pulsars^ 

PPTA 

Pai'kes 

64 

Australia 

20 


Effelsberg 

100 

Germany 



Lovell 

76.2 

UK 


EPTA 

Nancay 

94b 

France 

22 


Sardinia 

64 

Italy 



Westerbork 

96*’ 

Netherlands 


NANOGrav 

Ai'ecibo 

GET 

305 

100 

USA 

22 


Notes: ^ We present the number of pulsars that have been timed for more 
than 5 years for each project here (see table 3 in ref. 1221 for more informa¬ 
tion). It is worth pointing out that a number of pulsars were recently added 
to the timing arrays and some pulsars were dropped from these arrays. The 
total number of pulsars that are currently being observed for IPTA is 70, of 
which 13 ai‘e timed by two timing annys and 10 by all three airays. 

^ Values of circular-dish equivalent diameter ED 



Dec=-90 


Figure 1 Distribution of millisecond pulsars on the sky in equatorial co¬ 
ordinates, containing 70 IPTA pulsars (red stars) and all presently known 
pulsars that have pulse periods P < 15 ms and P < 10^*^ ss“* (blue dots; as 
found in the ATNF Pulsar Catalogue version 1.53 1291 ) 


In this paper, we first introduce the basics of pulsars and 
pulsar timing techniques in section |2] and [3 respectively. In 
section |4] we highlight some major noise sources that affect 
PTA’s sensitivities to GWs. In section |5] we describe how a 
PTA responds to GWs. Section |6] contains a review on GW 
sources and related results derived from the latest analysis of 
PTA data. Einally we discuss future prospects in sectionEl 

2 Pulsars 

Neutron stars are born as compact remnants of core collapse 
supernovae during the death of massive main sequence stars 
with masses of around 8-25 Mq. A pulsar is a highly mag¬ 
netized, spinning neutron star. It can be detected as it emits 
beams of electromagnetic radiation along its magnetic axis 
that is misaligned with its rotational axis. Such beams of radi¬ 
ation sweep over the Earth in the same way lighthouse beams 
sweep across an observer, leading to pulses of radiation re¬ 
ceived at the observatory. Because of their exceptional rota¬ 
tional stability, pulsars are powerful probes for a wide range 
of astrophysical phenomena. As mentioned above, long-term 
timing observations of PSR B1913-I-16 provided the first ob¬ 
servational evidence of the existence of GWs. Timing obser¬ 
vations of PSR B1257-1-12 led to the first confirmed discovery 
of planets outside our solar system ESI- The double pulsar 
system PSR J0737-3039A/B, with both neutron stars having 
been detected as radio pulsars II31II32L enabled very stringent 
tests of general relativity and alternative theories of gravity in 
the strong-field regime ll^ . 

The first pulsar was discovered by Bell and Hewish in 
1967 134). Since then over 2500 pulsars have been discov- 
erec0, with spin periods ranging from about 1 millisecond to 
10 seconds. It is generally believed that pulsars are born with 


*http://WWW.atn£.csiro.au/research/pulsar/ppta/ 

^ http: //WWW .epta.eu.org/ 

’http://nanograv.org/ 
http: //WWW . ipta4gw. org/ 

^See the ATNF Pulsar Catalogue website http://www. atn£.csiro.au/research/pulsar/psrcat/for up-to-date information. 
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periods of order tens of milliseconds but quickly (~ 10 Myr) 
spin down (because of the loss of rotational energy) to periods 
of order seconds. This energy loss could be due to a variety of 
mechanisms, such as the magnetic dipole radiation, emission 
of relativistic particle winds and even GWs. As pulsars spin 
down, they eventually reach a point where there is insufficient 
energy to power electromagnetic radiation. However, for pul¬ 
sars in binary systems, it is very likely that pulsars are spun 
up as mass and angular momentum are accreted from their 
stellar companions. Such an accretion process is observed in 
X-ray binaries. These pulsars, usually named as millisecond 
pulsars, have spin periods of about several milliseconds and 
much lower spin-down rates. Figure |2] shows the distribution 
of all known pulsars in the period-period derivative diagram. 



Figure 2 The period (P) vs. period derivative (P) diagram for all known 
pulsars (black dots; data taken from the ATNF Pulsai' Catalogue version 
1.53). Red(grey) circles represent the millisecond pulsars cun'ently being 
timed by the IPTA project. Also shown are lines of constant characteristic 
age T = P/2P (dash) assuming that pulsars are spinning down solely because 
of magnetic dipole radiation [33, and of constant infen'ed surface magnetic 
field Bq = 3.2 X 10^^ VpPGauss (solid; 1361 ). Two distinct populations are 
apparent in this diagram: (a) normal pulsars, with P ~ 0. seconds and 
Bq ~ 10^ ^-10^^ Gauss; (b) millisecond pulsars, with P ~ 3 milliseconds 
and P ~ 10“^*^ ss“^ 

3 Pulsar timing techniques 

Much of the science based on pulsar observations makes 
use of the “pulsar timing” technique which in¬ 

volves measurement and prediction of pulses’ times of arrival 
(TOAs). Individual pulses are generally not useful in this re¬ 
gard as they are unstable and mostly too weak to observe. The 
average pulse profile over a large number of pulses is stable 
for a particular pulsar at a given observing wavelength, and 
therefore very suitable for timing experiments. 

The first step in pulsar timing is to measure the topocen- 
tric pulse arrival times with clocks local to the radio obser¬ 
vatories. Data collected with the telescope are de-dispersed 


to correct for frequency-dependent dispersion delays due to 
the ionized interstellar medium. These data are then folded 
with the period derived from previous observations to form 
the mean pulse profile. This profile is then correlated with a 
standard template, either an analytic function or simply a very 
high signal-to-noise ratio observation, to record the pulse ar¬ 
rival time at the observatory. 

The measured TOAs are further transformed to the pulse 
emission time via a timing model, from which the pulse phase 
of emission is computed. The rotational phase (j){t) of the pul¬ 
sar as a function of time t (measured in an inertial reference 
frame of the pulsar) can be represented as a Taylor series: 

0(0 = 0(fo) + f{t - to) + -f{t - to)^ + ■ ■■ , (1) 

where to is an arbitrary reference time, / = d0/df is the 
spin frequency, and / is the frequency derivative. A number 
of corrections are applied when converting the topocentric 
TOAs to the pulsar frame. Such corrections include: 

1. Clock corrections, which account for differences in the 
observatory time and a realization of Terrestrial Time 
(e.g., the International Atomic Time). 

2. pulse delay induced by Earth’s troposphere. 

3. the Einstein delay, i.e., the time dilation due to changes 
in the gravitational potential of the Earth, the Earth’s 
motion, and the secular motion of the pulsar or that of 
its binary system. 

4. the Roemer delay, i.e., the vacuum light travel time be¬ 
tween the observatory and the solar system barycenter, 
and for pulsars in binaries between the pulsar and the 
binary system’s barycenter. 

5. Shapiro delays, i.e., gravitational time delays due to 
the solar system objects and if applicable the pulsar’s 
companion. 

6. Dispersion delays caused by the interstellar medium, 
the interplanetary medium and the Earth’s ionosphere. 

The timing model, which describes the above corrections 
and the pulsar’s intrinsic rotational behavior, predicts the ro¬ 
tational phase of the pulsar at any given time as observed 
from the solar system barycenter. Basic parameters of a pul¬ 
sar timing model include the spin period, spin-down rate, 
right ascension and declination of the pulsar, the dispersion 
measure (discussed later in section Si, and Keplerian orbital 
parameters if the pulsar is in a binary system. Measured 
TOAs are compared with predictions based on the timing 
model, and the differences are called timing residuals. The 
(pre-fit) timing residual for the i-th observation is calculated 
as 139): 
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where Ni is the nearest integetQ to each 0(f,). One can see 
the key point in pulsar timing is that every single rotation of 
the pulsar is unambiguously accounted for over long periods 
(years to decades) of time. 

A linear least-squares fitting procedure is carried out to 
obtain estimates of timing parameters, their uncertainties and 
the post-fit timing residuals. In practice, this is done itera¬ 
tively: one starts from a small set of data and only includes 
the most basic parameters (with values derived from previous 
observations) so that it is easier to coherently track the rota¬ 
tional phase. Parameter estimates are then improved by mini¬ 
mizing the timing residuals and additional parameters can be 
included for a longer data set. 

The fitting to a timing model and analysis of the tim¬ 
ing residuals can be performed with the pulsar timing soft¬ 
ware package TEMP02 ||32ll39]l^, which is freely avail¬ 
able on the internet for downloaqj. More recently, an alter¬ 
native method based on Bayesian inference was also devel¬ 
oped 041II42I . 

4 Noise sources in pulsar timing data 

Timing residuals generally come from two groups of con¬ 
tribution: (1) un-modelled deterministic processes, e.g., an 
unknown binary companion or a single-source GW; and (2) 
stochastic processes, e.g., the intrinsic pulsar spin noise and 
a GW background. Before the discussion of GW detection 
with PTAs, we briefly discuss some major noise processes in 
pulsar timing data here. 

Radiometer noise 

Radiometer noise arises from the observing system and the 
radio sky background (including the atmosphere, the cos¬ 
mic microwave background and synchrotron emission in the 
Galactic plane). It can be quantified as llTSl : 

w / W 

‘^”‘' ~5/A~5™eanV2A^ Vp-W’ 

where W and P are the pulse width and period respectively, 

5 IN is the profile signal-to-noise ratio, 5sys is the system 
equivalent flux density which depends on the system tem¬ 
perature and the telescope’s effective collecting area, 5 mean 
the pulsar’s flux density averaged over its pulse period, Av 
and tint are the observation bandwidth and integration time 
respectively. Radiometer noise can be reduced by using low- 
noise receivers, observing with larger telescopes, and increas¬ 
ing observing time and bandwidth. One reason to fold indi¬ 
vidual pulses to obtain an average profile is to reduce the ra¬ 
diometer noise; the reduction is equal to the square root of 
the number of pulses folded. Radiometer noise is an additive 
Gaussian white noise and is formally responsible for TOA 
uncertainties. From equation (O, one can see that bright, fast 
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spinning pulsars with narrow pulse profiles allow the highest 
timing precision. 

Pulse jitter noise 

Pulse jitter manifests as the variability in the shape and ar¬ 
rival phase of individual pulses. The mean pulse profile is 
an average over a large number of single pulses. Although 
it is stable for most practical purposes, there always exists 
some degree of stochasticity in the phase and amplitude of 
the average pulse profile. Pulse jitter noise is intrinsic to the 
pulsar itself, and thus can only be reduced by increasing ob¬ 
serving time, i.e., averaging over more single pulses. Pulse 
jitter is also a source of white noise, which has been found 
to be a limiting factor of the timing precision for a few very 
bright pulsars Il43ll44l . For future telescopes such as FAST 
and SKA, jitter noise may dominate over the radiometer noise 
for many millisecond pulsars 1261 . Improvement in the tim¬ 
ing precision for the brightest pulsar PSR J0437-4715 was 
recently demonstrated with the use of some mitigation meth¬ 
ods for pulse jitter noise 145II46L 

“Timing noise” 

It has long been realized that timing residuals of many pul¬ 
sars show structures that are inconsistent with TOA uncer¬ 
tainties EHH. Such structures are collectively referred to 
as timing noise. Timing noise is very commonly seen in nor¬ 
mal pulsars and has also become more prominent in a number 
of millisecond pulsars as timing precision increases and the 
data span grows. The exact astrophysical origins of timing 
noise are not well understood. It is mostly suggested to be 
related to rotational irregularities of the pulsar and therefore 
it is also usually called as spin noise II491I50II . For millisec¬ 
ond pulsars, power spectra of timing residuals can typically 
be modelled as the sum of white noise and red noise. For red 
noise, a power law spectrum with a low-frequency turnover 
appears to be a good approximation (see, e.g., figure 11 in 
ref. IfTTlI for analyses of 20 PPTA pulsars). 

The presence of red noise results in problems to the pul¬ 
sar timing analysis as the standard least-squares fitting of a 
timing model assumes time-independent TOA errors. Bland- 
ford et al. Il47ll analytically showed the effects of timing noise 
on the estimates of timing model parameters and suggested 
the use of the noise covariance matrix to pre-whiten the data 
for improved parameter estimation. More recently. Coles et 
al. ED developed a whitening method that uses the Cholesky 
decomposition of the covariance matrix of timing residuals to 
whiten both the residuals and the timing model. By doing so, 
noise in the whitened residuals is statistically white and the 
ordinary least-squares solution of a timing model can be ob¬ 
tained. van Haasteren & Levin developed a Bayesian 


®Here note that 0(t) is measured in turns equal to 27t radians. 
2 www.sf.nel/projects/tempo2/ 
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framework that is capable of simultaneously estimating tim¬ 
ing model parameters and timing noise spectra. 

Dispersion measure variations 

Because of dispersion due to the interstellar plasma, pulses at 
low frequencies arrive later than at high frequencies. Specif¬ 
ically, this dispersion delay is given by: 

ADM^(4.15ms)DMvc2^^, (4) 

where vghz is the radio frequency measured in GHz, and dis¬ 
persion measure (DM, measured in pc cm“^) is the integrated 
column density of free electrons between an observer and 
a pulsar. Because of the motion of the Earth and the pul¬ 
sar relative to the interstellar medium, the DM of a pulsar is 
not a constant in time. Such DM variations introduce time- 
correlated noise in pulsar timing data. 

Pulsars in a timing array are usually observed quasi- 
simultaneously at two or more different frequencies. For ex¬ 
ample, the PPTA team observes pulsars at three frequency 
bands - 50 cm (-700 MHz), 20 cm (-1400 MHz), and 10 
cm (-3100 MHz) - during each observing session (typically 
2-3 days). This makes it possible to account for DM varia¬ 
tions and thus reduce the associated noise with various meth¬ 
ods lEMa. The method currently being used by the PPTA 
is described in Keith et al. JSl, which was built on a previ¬ 
ous work by You et al. 15^ . In this method timing residu¬ 
als are modelled as the combination of a (radio-)wavelength- 
independent (i.e., common-mode) delay and the dispersion 
delay. Both components are represented as piece-wise lin¬ 
ear functions and can be estimated through a standard least- 
squares fit. A key feature of this method is that GW signals 
are preserved in the common-mode component. Ultra wide¬ 
band receivers are in development by various collaborations 
in order to better correct for noise induced by DM variations 
(along with other benefits such as increasing timing precision, 
studies of interstellar medium, and etc.). 

Interstellar scintillation 

Interstellar scintillation refers to strong scattering of radio 
waves due to the spatial inhomogeneities in the ionized in¬ 
terstellar medium Et), analogous to twinkling of stars due to 
scattering in the Earth’s atmosphere. There are multiple ef¬ 
fects associated with interstellar scintillation that cause time- 
varying delays in measured TOAs, with the dominant one be¬ 
ing pulse broadening from multipath scattering 058II59I . Var¬ 
ious mitigation techniques have been developed for this type 
of noise EMI]. Generally speaking, noise induced by in¬ 
terstellar scintillation is a Gaussian white noise, and can be 
reduced by increasing observing time and bandwidth. For 
millisecond pulsars that are observed at current radio fre¬ 
quencies by PTAs, the effects of scintillation are predicted to 
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be small. However, when pulsars are observed at lower fre¬ 
quencies, or more distant (and more scintillated) pulsars are 
observed, these effects can become more important II59II62I . 

Correlated noise among different pulsars 

The above noise processes are generally thought to be uncor¬ 
related among different pulsars. In a PTA data set that we 
hope to detect GWs, some correlated noise may be present. 
For example, (1) instabilities in Terrestrial Time standards 
affect TOA measurements of all pulsars in exactly the same 
way, i.e., clock errors result in a monopole signature in a PTA 
data set; (2) the solar system ephemerides, which provide ac¬ 
curate predictions of the masses and positions of all the major 
solar system objects as a function of time, are used to convert 
pulse arrival times at the observatory to TOAs referenced at 
the solar system barycenter. Imperfections in the solar sys¬ 
tem ephemerides induce a dipole correlation in a PTA data 
set. Indeed, it has been demonstrated that PTAs can be used 
(1) to search for irregularities in the time standard and thus 
to establish a pulsar-based timescale ll6Tl . and (2) to measure 
the mass of solar system planets ||6^ . 

5 Gravitational waves and pulsar timing ar¬ 
rays 

The effects of GWs in a single-pulsar data may be indistin¬ 
guishable from those due to noise processes as discussed in 
the previous subsection. Indeed, even without any such noise, 
GWs that have the same features as those due to uncertainties 
in the timing model parameters would still be very difficult 
to detect with only one pulsar. Therefore analysis of single¬ 
pulsar data can be best used to constrain the strength of po¬ 
tential GWs Il^l66l. 

A PTA is a Galactic-scale GW detector. If one wishes to 
have an analogy to a laser interferometer, pulsars in the tim¬ 
ing array are “test masses”; pulses of radio waves act as the 
laser; and the pulsar-Earth baseline is a single “arm”. Mil¬ 
lisecond pulsars in our Galaxy, typically ~kpc (or thousands 
of light years) away, emit radio waves that are received at the 
telescope with extraordinary stability. A GW passing across 
the pulsar-Earth baseline perturbs the local spacetime along 
the path of radio wave propagation, leading to an apparent 
redshift in the pulse frequency that is proportional to the GW 
strain amplitude. Fet us first consider the special case where 
a linearly polarized GW propagates in a direction perpendic¬ 
ular to the pulsar-Earth baseline, the resulting timing residual 
is given by: 

r{t) =1 h(t -+ T)dT, (5) 

Jo c 

where L is the pulsar distance and we adopt the plane wave 
approximations With the definition of cL4(f)/df = h(t), the 


^For sources that are close enough (< 100 Mpc), it may be necessaiy to consider the curvature of the gravitational wavefront. This, in principle, would 
allow luminosity distances to GW sources to be measured via a pai'allax effect (67). 
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timing residual takes the following form: 

r{t) = AA(f) = A{t) - Ait - -). (6) 

c 

Here we can see that A(f) results from the GW induced space- 
time perturbation incident on the Earth (i.e., the Earth term), 
and Ait- depends on the GW strain at the time of the radio 
wave emission (i.e., the pulsar term). Typical PTA observa¬ 
tions have a sampling interval of weeks and span over ~ 10 yr, 
implying a sensitive frequency range of ~ 1-1 00 nHz. There¬ 
fore PTAs are sensitive to GWs with wavelengths of several 
light years, much smaller than the pulsar-Earth distance. 

In the general case where a GW originates from a direction 
Q, the induced timing residuals can be written as: 

rit, Q) = E+(Q)AA^(r) + Ex(Q)AAx(f), (7) 

where E+(Q) and F^iCl) are antenna pattern functions as 
given by Il68l : 

^ 77r~ -^{(1+sin^^)cos2b,,cos[2(a-Q'p)] 

4( 1 - cos 0) 

- sin 2(5 sin 26p cos(a; - ap) + cos^ 5(2 - 3 cos^ 5p)) (8) 

Ex(Q) = ——-^{cos 6 sin 26p sin((3' - ap) 

Z\i — cos u) 

-sin5cos^5pSin[2(a-Q-p)]), (9) 

where cos 6 - cos 6 cos 6p cos(q: - ap) + sin 6 sin 6p, 9 is the 
opening angle between the GW source and pulsar with re¬ 
spect to the observer, and a iap) and 6 i6p) are the right as¬ 
cension and declination of the GW source (pulsar) respec¬ 
tively. The source-dependent functions AA+_x(0 in equation 
© are given by: 

AA+,x(f)=A+,x(r)-A+,x(fp) (10) 

f„ = f - (1 - COS0) — . (11) 

c 

The forms of A+(f) and Ax(f) depend on the type of source 
that we are looking for. 
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larger than the original result of ref. Q. This is because ^i9ij) 
is normalized to 1 for the autocorrelation of the stochastic 
background induced timing residuals for a single pulsar. In 
Eigure|3]the correlation function takes a value of 0.5 at zero 
angular separation as the autocorrelation due to pulsar terms 
is neglected. 



Figure 3 The Hellings-Downs cuiwe, which depicts the expected corre¬ 
lation in timing residuals due to an isotropic stochastic background, as a 
function of the angular separation between pairs of pulsars. 

6 Gravitational wave sources and recent obser¬ 
vational results 

Potential signals that could be detectable for PTAs include: 

(1) stochastic backgrounds. The primary target is that formed 
by the combined emission from numerous binary supermas- 
sive black holes distributed throughout the Universe. Back¬ 
ground signals from cosmic strings (e.g., 16911 ) and inflation 
(e.g., 070117 II ) have also been studied in the context of PTAs; 

(2) continuous waves, which can be produced by individual 
nearby binaries; (3) bursts with memory associated with bi¬ 
nary black hole mergers; (4) and all other GW bursts. Be¬ 
low we give a brief overview of these sources and summarize 
some recent astrophysical results. 


The Hellings-Downs curve 

An isotropic stochastic background will produce a correlated 
signal in PTA data sets. Such a correlation uniquely depends 
on the angular separation between pairs of pulsars, as given 
by©: 


mj) 


3(1- cos 9ii) [(1 - cos 6ij) 

2 - 2 -- 2 - - 

l(l-cos%) il+6ij) 

4 2 ^ 2 ’ 


( 12 ) 


where Ojj is the angle between pulsars i and j, and 5,^ is 1 for 
i - j and 0 otherwise. Eigure [3 shows the famous Hellings- 
Downs curve as given by equation (fT2l i - it is a factor of 3/2 


Stochastic backgrounds 

The stochastic background from the cosmic population of su- 
permassive binary black holes has been the most popular tar¬ 
get for PTA efforts. Generally speaking the signal amplitude 
depends on how frequently these binaries merge in cosmic 
history and how massive they are. Both of these quantities 
are poorly constrained observationally. Assuming that all bi¬ 
naries are in circular orbits and evolve through gravitational 
radiation only, the characteristic amplitude spectrum of this 
background is given by 11721 - 1781 : 


, ( 13 ) 
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where Ayr is the dimensionless amplitude at a reference fre¬ 
quency /yr = 1 yr The fraction of the cosmological critical 
energy density (per logarithmic frequency interval) contained 
in the GW background is related to the amplitude spectrum 
through Qgw(/) = where Hq is 

the Hubble constant. Various models predict a similar range 
of Ayr, most likely to be ~ 10“'^ (see, e.g., ElHnl). An 
exception is the recent model in ref. llTSll whose prediction 
is two to five times higher. Recent studies that include the 
effects of environmental coupling and orbital eccentricities 
indicate a reduced signal at below ~ 10 nHz ITTOHSTII . 

Jenet et al. Il82l suggested that it is possible to make a de¬ 
tection of the binary black hole background if 20 pulsars are 
timed with a precision of ~ 100 ns over > 5 years. Three 
PTAs have searched for such a background signal assuming 
it is isotropic, leading to increasingly more stringent upper 
limits on the background strength 053ll83fl87l . The most con¬ 
straining limit published to date (Ayr < 2.4 x 10“'^) comes 
from the PPTA collaboration by Shannon et al. IISTI . As 
shown in Figure |4] this limit ruled out the most optimistic 
model of iTTSl with 90% confidence and is in tension with 
other models at ~ 50% confidence. 


Recently methods have also been proposed to search for a 
more general anisotropic background signal ||9TI494|| . Using 
the 2015 EPTA data, Taylor et al. ||95]| placed constraints on 
the angular power spectrum of the background from circular, 
GW-driven supermassive black hole binaries and found that 
the data could not update the prior knowledge on the angular 
distribution of a GW background. 


Continuous waves 

Individual supermassive binary black holes, especially the 
most nearby and/or massive ones, could provide good oppor¬ 
tunities for detection of continuous waves. Unlike compact 
binaries in the audio band, supermassive binary black holes 
detectable for PTAs are mostly in the early stage of inspiral 
and therefore emit quasi-monochromatic waves. For an in- 
spiralling circular binary of component masses nii and m 2 , 
the GW strain amplitude is given by ||9^ : 

ho = 2- 2 ---, (14) 

di 



Figure 4 Upper limits on the fractional energy density of the GW back¬ 
ground O.Qyy{f), at a frequency of 2.8 nHz, as compared against various 
models of the supermassive binary black hole background. The red(grey) 
solid and dashed lines show the probabilities Pr(nGw) that a GW background 
with energy density HowC/) exists given the PPTA data, assuming Gaus¬ 
sian and non-Gaussian statistics respectively. Three vertical lines mark the 
95% confidence upper limits published by NANOGrav 1531 . EPTA 1851 and 
PPTA (87), where the times next to these lines correspond approximately 
to the observing spans of the data sets. The shaded region is ruled out with 
95% confidence by the PPTA data. Gaussian cuiwes represent the probability 
density functions pM(f^Gw) given different models: a merger-driven model 
for growth of massive black holes in galaxies GIl, a synthesis of empirical 
models ESI, semi-analytic models (SMBH model bEIl) based on the Mil¬ 
lenium dark matter simulations I88II891 and a distinct model for black hole 
growth (SMBH model 2; 1901 ). Figure taken from ref. l87l . 


where di is the luminosity distance of the source, and Me is 
the chirp mass defined as Me = Mrf^^, with M - mi + m 2 
the total mass and rj - mim 2 fM^ the symmetric mass ratio. 
After averaging over the antenna pattern functions given by 
equations (I8]l9]l and the binary orbital inclination angle, the 
Earth-term timing residuals induced by a circular binary is 
~ hol2nf. 

Before the establishment of major PTAs, Jenet et al. ll^ 
developed a framework in which pulsar timing observations 
can be used to constrain properties of supermassive binary 
black holes and applied the method to effectively rule out the 
claimed binary black hole system in 3C 66B ||971 . In recent 
years growing efforts have gone into investigating the detec¬ 
tion prospects Il98l4102l of, and designing data analysis meth¬ 
ods 11103111091 for, continuous waves. 

Yardley et al. Ill 101 calculated the first sensitivity curve 
of a PTA to this type of sources using an earlier PPTA data 
set presented in ref. Ill 1 II . Recently both PPTA II1081 and 
NANOGrav 111 121 conducted searches for continuous waves 
in their corresponding real data sets. Because of its excellent 
data quality, PPTA has achieved by far the best sensitivity for 
continuous waves. Pigure|5] shows a sky map of PPTA’s sen¬ 
sitivities to circular binary black holes of chirp mass lO^M© 
and orbital period of ~6 years 11081 . Unfortunately most 
nearby galaxy clusters or binary black hole candidates are 
within the less sensitive sky region. Zhu et al. nm also 
presented very stringent upper limil^ on the GW strain am¬ 
plitude (ho < 1.7 X 10 at 10 nHz) and on the local merger 
rate density of supermassive binary black holes (fewer than 
4 X lO^^Mpc^^Gyr^* forM^ > IO'^Mq). 


More recently, the EPTA group obtained slightly better limits for / < 10 nHz and comparable results at higher frequencies Eni- 
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Dec=90 



Figure 5 Sky distribution of luminosity distance (in Mpc) out to which a 
circular binary black hole of chirp mass 10^M© and orbital period of 5 nHz 
could be detected with the PPTA data set. Sky locations of the 20 PPTA 
pulsars are indicated by White diamonds mark the location of possible 
supermassive black hole binary candidates or nearby clusters: Virgo (16.5 
Mpc), Fornax (19 Mpc), Norma (67.8 Mpc), Perseus (73.6 Mpc), 3C 66B 
(92 Mpc), Coma (102 Mpc) and OJ287 (1.07 Gpc). Figure originally pub¬ 
lished in ref. (TOD. 

Gravitational wave memory 

A GW memory is a permanent distortion in the spacetime 
metric 111 14111151 . Such effects can be produced during merg¬ 
ers of supermassive binary black holes and cause instanta¬ 
neous jumps of pulse frequency. For a single pulsar, this 
is indistinguishable from a glitch event. With a timing ar¬ 
ray, GW memory effects can be searched for as simultaneous 
pulse frequency jumps in all pulsars when the burst reaches 
the Earth. It has been suggested that GW memory signals 
are in principle detectable with current PTAs for black hole 
mass of 10^ Mq within a redshift of 0.1 0116411191 . How¬ 
ever, the event rate is highly uncertain. Current estimates 
are very pessimistic, predicting only 0.03 to 0.2 detectable 
events every 10 years for future PTA observations based on 
the SKA 0101111201 . Actual searches in existing PTA data 
sets - see ref. 01211 for PPTA and ref. 01220 for NANOGrav 
- have set upper limits on the memory event rate, which re¬ 
main orders of magnitude above theoretical expectations. 

Gravitational wave bursts 

Potential burst sources of interest to PTAs include the forma¬ 
tion or coalescence of supermassive black holes, the periap- 
sis passage of compact objects in highly elliptic or unbound 
orbits around a supermassive black hole II123L cosmic (su- 
perjstring cusps and kinks 11124141261 . and triplets of super¬ 
massive black holes M127II . Finn (& Lommen M123II developed 
a Bayesian framework for detecting and characterizing burst 
GWs (see also ref. 111281 ). Zhu et al. II1091 recently proposed 
a general coherent (frequentist) method that can be used to 
search for GW bursts with PTAs. No specific predictions 
have been made for detecting GW bursts with PTAs in the 
literature. There has been no published results on searches 
for bursts using real PTA data. 


7 Summary and future prospects 

With their existence predicted by general relativity one cen¬ 
tury ago, GWs have not yet been directly detected. However, 
it is widely believed that we are on the threshold of open¬ 
ing the gravitational window into the Universe. In the au¬ 
dio band (from 10 to several kilo Hz), 2nd-generation laser 
interferometers such as Advanced LIGO are about to start 
scientific observations and a detection of signals from com¬ 
pact binary coalescences (e.g., binary neutron star inspirals) 
is likely within a few years (see the article by Reitze et al. 
in this issue). In the nanohertz frequency range, PTA exper¬ 
iments have achieved unprecedented sensitivities and started 
to put serious constraints on the cosmic population of super¬ 
massive black hole binaries. 

The sensitivity of a PTA can be improved by a) increas¬ 
ing the data span and observing cadence, b) including more 
pulsars in the array, and c) reducing the noise present in the 
data. Regarding the first two factors, the combination of data 
sets from three PTAs to a single IPTA data set offers the 
most straightforward benefit. Other ongoing efforts include 
optimization of observing strategies, searches for millisec¬ 
ond pulsars, characterization of various noise processes and 
corresponding mitigation methods, and development of ad¬ 
vanced instrumentations. In the longer term, pulsar timing 
observations with FAST and SKA will provide advances in 
all aspects of PTA science, not only leading to the detection 
of GWs but also allowing detailed studies of the nanohertz 
gravitational Universe. 
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